Non-affine response: jammed packings versus spring networks 
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Abstract. - We compare the elastic response of spring networks whose contact geometry is derived 
from real packings of frictionless discs, to networks obtained by randomly cutting bonds in a highly 
connected network derived from a well-compressed packing. We find that the shear response of 
packing-derived networks, and both the shear and compression response of randomly cut networks, 
are all similar: the elastic moduli vanish linearly near jamming, and distributions characterizing 
the local geometry of the response scale with distance to jamming. Compression of packing-derived 
networks is exceptional: the elastic modulus remains constant and the geometrical distributions 
do not exhibit simple scaling. We conclude that the compression response of jammed packings is 
anomalous, rather than the shear response. 



The jamming transition governs the onset of rigidity 
in disordered media as diverse as foams, colloidal suspen- 
sions, granular media and glasses [1]. While jamming in 
general is controlled by a combination of density, shear 
stress and temperature, most progress has been made for 
frictionless soft spheres that interact through purely re- 
pulsive contact forces, and that are at zero temperature 
and zero load [2-7]. This simple model applies to static 
foams or emulsions [8,9], and represents a simplified ver- 
sion of granular media, if one ignores friction [10,11] and 
nontrivial grain shapes [12-15]. 

From a theoretical point of view, this model is ideal 
for two reasons. First, it exhibits a well defined jam- 
ming point, "point J", which in the limit of large system 
sizes, occurs at a well-defined density <j) = <f) c [2]. Here the 
system is a disordered packing of frictionless undeformed 
spheres, which is marginally stable and isostatic, i.e., its 
contact number (average number of contacts per particle) 
z equals Zi so = 2c? in d dimensions [2,16]. Second, in re- 
cent years it has been uncovered that the mechanical and 
geometric properties of such jammed packings exhibit a 
number of non-trivial power law scalings as a function of 
the distance to the jamming point: (1) The excess con- 
tact number Az :— z — z- lso scales as (4>— fic) 1 ^ 2 [2,6,9,10]; 
(2) The ratio of shear (G) and bulk (K) elastic moduli 
vanishes at point J as G/K ~ Az [2]. 

The latter behavior — a shear rigidity which becomes 



much smaller than the compression modulus as the jam- 
ming point is approached - is in many ways surprising. 
It also differs markedly from what is found in two sim- 
plified models of jammed systems, effective medium the- 
ory (EMT) and random clastic networks, as is illustrated 
schematically in fig. 1 for the simple case of harmonic 
particles. EMT predicts that the elastic moduli vary 
smoothly through the isostatic point where Az = and 
that the moduli are of order of the local spring constant 
k. This is because effective medium theory is essentially 
"blind" to local packing considerations and isostaticity. 
Thus, besides failing to capture the vanishing of G near 
jamming, its prediction for the bulk modulus fails spectac- 
ularly as well: it predicts finite rigidity below isostaticity. 

The failure of EMT to describe elasticity near jamming 
motivated earlier suggestions that elasticity of jammed 
packings might be captured by random networks of springs 
— this problem is known as rigidity percolation [8,17-19]. 
However, in such random spring networks, both G and K 
are expected to go to zero as kAz, as fig. lc illustrates [17]. 

Thus, while from the point of view of effective medium 
theory the shear rigidity of jammed packings behaves 
anomalously, from the point of view of rigidity percolation, 
the compression modulus behaves unexpectedly What 
sets jammed packings apart from either of these two lim- 
iting models? How to understand the difference in terms 
of the local packing or response? Is the difference with 
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Fig. 1: Schematic comparison of the variation of shear (G) and 
bulk (K) elastic moduli as function of distance to jamming, 
(a) In effective medium theory, all elastic moduli are simply 
of the order of the local spring constant k, and moreover, the 
theory does not account for whether the packing is rigid or 
not. (b) In jammed packings of harmonic particles, the bulk 
modulus K remains constant down to the jamming transition, 
where it vanishes discontinuously, whereas the shear modulus 
G vanishes linearly in Az. (c) In random networks of elastic 
springs, both elastic moduli vanish linearly with Az. 



rigidity percolation visible in the scaling behavior of the 
response of packings? These are issues we aim to clarify 
in this paper. Our approach will hinge on characterizing 
the elastic response at the level of the bonds. After all, 
the elastic moduli characterize changes in elastic energy 
AE under deformations, and AE simply is a sum of the 
changes in elastic energy of all contacts (bonds) in the 
system. 

By probing the nature of the local response of packing- 
derived and randomly cut networks, we find that we can 
distinguish two cases. In the "generic" case, all geometri- 
cal characterizations exhibit simple scaling and the elastic 
moduli scale as Az — this describes shear and bulk de- 
formations of randomly cut networks, as well as shear de- 
formations of packing-derived networks. Packing derived 
networks under compression form the "exceptional" case: 
the fact that the compression modulus remains of order 
k near jamming is reflected in the fact that various char- 
acteristics of the local displacements do not exhibit pure 
scaling. We connect these findings to recent theoretical 
work by Wyart [20,21]. 

Linear Response. — All numerical results presented 
in this paper concern quasistatic linear response of sys- 
tems to global shear or comprcssional forcing. First we 
generate, for a range of pressures, ensembles of 50 two- 
dimensional jammed packings of 1024 frictionless particles 
with one-sided harmonic forces (k = 1) using a Molecu- 
lar Dynamics simulation (for details, see [22]). Our linear 
response calculations are based on the dynamical matrix. 
We decompose, for linear deformations, the relative dis- 
placement Uij of neighboring particles i and j in compo- 
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Fig. 2: Two families of spring networks — see text for details. 



Tij connects the centers of particles i and j. In these terms 
the change in energy takes a simple form [5,23], 



AE = 



2^2 



fi. 



(1) 



The dynamical matrix Mij, a p is obtained by rewriting 
eq. (1) in terms of the independent variables, Ui t0t , as 



AE = ljM t j, a p u i>a Uj.p . 



(2) 



Here are contact forces, k denotes the stiffness of the 
springs 1 , M. is a dN x dN matrix with N the number 
of particles, indices a, (3 label the coordinate axes, and 
the summation convention is used. The dynamical matrix 
contains all information on the elastic properties of the 
system, and in particular describes the linear response to 
external forcing as [7, 24] : 



(3) 



Two Families of Spring Networks. — We start by 
noting that the analysis of the linear response of jammed 
packings of particles with one-sided harmonic interactions 
(fig. 2a) is exactly equivalent to that of networks of appro- 
priately loaded harmonic springs (fig. 2b), with the nodes 
of the network given by the particle centers and the geom- 
etry and forces of the spring network determined by the 
force network of the packing. 

In all that follows, we ignore the pre-stress term 
j.- which is subdominant near jamming — we have 



fi. 



nents parallel 



1 For our harmonic potential, k = 1 for each contact, but the 
procedure works equally well for more general potentials, for which 
is simply the value of the second derivative of the potential, 
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and perpendicular (ltj_) to ry, where evaluated at the initial distance n 
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Fig. 3: Shear (diamonds) and bulk (squares) elastic moduli 
as function of distance to jamming for (a) jammed harmonic 
packings (b) randomly cut spring networks. Both panels have 
the same linear fit for the shear modulus (blue line). Inset: 
enlargement of the behavior of K and G near z = 4. 



checked that its exclusion does not affect the results [6]. 
The system without pre-stress is equivalent to a "neutral" 
spring network where all contacts are replaced by springs 
at their equilibrium length so that /y = for all contacts 
(fig. 2c). For such a neutral spring network the dynamical 
matrix becomes particularly simple, as its only non-zero 
elements are simply given by geometry and by the bond 
strengths k of each bond. 

We follow two routes to approach the (un)jamming tran- 
sition by lowering the contact number in the neutral net- 
works. In the first route, we map jammed packings un- 
der increasingly low pressure (fig. 2a) to neutral, packing- 
derived spring networks (fig. 2c) — the geometry of these 
networks and the contacts that are broken when point J 
is approached are thus set by real packings that were cre- 
ated using the MD protocol described in Rcf. [22]. In 
the second route, we start from a neutral spring network 
that is obtained from a heavily compressed jammed pack- 
ing for which z ss 5.98. We then create randomly cut 
networks with lower contact number by randomly remov- 
ing springs [25], making sure that we do not create local 
disconnected patches and that each node in the network 
remains connected by at least three springs (fig. 2d) — the 
geometry of these networks becomes increasingly random. 
Note that no relaxation is needed after removing springs 
because the neutral network has — in each contact. 

Elastic Moduli. — To analyze the linear response, 
we impose an infinitesimal strain deformation of order e, 
implemented by the appropriate changes in rest lengths of 
all bonds that cross the boundary of the periodic box of 
size LxL. This amounts to replacing Uy in eq. (1) by u.y — 
u]y , where denotes the appropriate shift of magnitude 
sL at bonds ij that cross the boundary, and is zero for 
interior bonds. Keeping track of this substitution in going 
from cq. (1) to eq. (3), these constant terms are taken to 
the right hand side, and thus act like an effective / ext [7] 
that is proportional to e. The response of the system to 
this shape or volume change of the box is then calculated 
by solving equation (3) for this effective external force. 

The moduli are extracted from the energy (eq. (2)) ac- 



K,G 



AE 
2V^ 



(4) 



for a uniform strain, e xx = e yy = e for compression, and 
£xy = £ for shear. Here V is the volume of the system. 

In fig. 3 we show the scaling of the elastic moduli G and 
K thus obtained, as a function of the contact number z for 
both packing-derived and randomly cut spring networks. 
The main point is that these two families exhibit different 
scaling behavior: for randomly cut networks, both mod- 
uli vanish as Az, while for the packing-derived networks 
only the shear modulus G goes to zero — the compression 
modulus K remains finite 2 . The behavior of the randomly 
cut networks is consistent with what is expected for rigid- 
ity percolation in random spring networks [17,18], while 
the behavior for packing-derived networks is in agreement 
with earlier data for jammed packings [2,7]. Hence, from 
the point of view of rigidity percolation, the anomaly in 
jammed packings is thus that the compression modulus 
K/k stays finite as Az — > 0. 

Note that the dynamical matrix of both types of net- 
works contains only geometric information about the net- 
work, since the spring constant k = 1 for each bond. Hence 
the crucial difference between packing-derived networks 
and randomly cut networks that is causing the bulk mod- 
ulus to be different must have a purely geometric origin. 

Nonaffinity of Response. — We will now connect 
the scaling of the elastic moduli to the local deformation 
field. One tool that we use to probe the degree of non- 
affinity of the response near point J is P{a), the proba- 
bility density function (PDF) of the displacement angles 
otij [7]. Here a denotes the angle between u»j and r^, or, 



tan on 



(5) 



In EMT, the displacements of the particles are pre- 
scribed by an affine deformation field. Affine compression 
corresponds to a uniform shrinking of the bond vectors, 



i.e. u± 



while 



-STij < 0: the corresponding 



P{ol) exhibits thus a delta peak at a — n. The effect of an 
affine shear on a bond vector depends on its orientation, 
and for isotropic random packings P(a) is flat. 

The results for packing-derived networks are shown in 
fig. 4ac. Note that far away from jamming, the PDFs are 
similar to the EMT predictions: a peak at tt under com- 
pression, and a flat PDF under shear. When approaching 
the unjamming transition, a peak at a = n/2 develops, 
which signifies that an increasingly large fraction of con- 
tacting particles mostly slide past each other. However, 
under shear, this peak is much more pronounced than un- 
der compression, and under compression the PDF retains 
a significant shoulder between 7r/2 and tt. 



2 Here and in what follows, Az = z — z c « z — Zi BO , where for 
packing-derived (randomly cut) networks z c = 4 (4.045). The dis- 
crepancy between z c and Zj so (see inset of fig. 3) for the randomly 
cut networks is not a finite size effect, but can be attributed to the 
precise cutting protocol. 
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Fig. 4: The PDF of the displacement angles P(a) for 
compression (a,b) and shear (c,d). The seven curves de- 
note, in order of decreasing peak height at a = n/2, 
z = 4.008, 4.027, 4.063, 4.14, 4.28, 4.74, 5.27. (a,c) For packing- 
derived networks, P(a) for compression and shear appear 
rather different. (b,d) For randomly cut networks, P(a) de- 
velops the same peaked structure when z ^ 4 for compression 
and shear. Insets: The width of the peaks (defined as the 
width of the interval containing the central 50% of the data), 
as a function of Az. The dotted lines indicate w ~ Az for all 
cases except compression of bead packings. 



The results from the randomly cut networks are shown 
in fig. 4bd: a strong peak develops in P(a) as Az de- 
creases, both for the response to shear and to compression. 
The relative displacements of contacting particles in re- 
sponse to compression thus signal an important difference 
between packing-derived networks and random networks. 

Scaling Arguments for Non Affinity. Wyart 
and coworkers have given arguments for estimating the en- 
ergies and local deformations of soft (low energy) modes 
starting from purely floppy (zero energy) modes [5,25]. 
They construct trial soft modes that are basically floppy 
modes, obtained by cutting bonds around a patch of size 
£* , and then modulating these with a sine function of 
wavelength £* to make the displacements vanish at the 
locations of the cut bonds. Here £* ~ 1/Az is a charac- 
teristic length scale [5-7,20]. In particular, for the local 
deformations (see fig. 5), they find [25] 



1 



3_ 



Az , 



(6) 



where symbols without indices ij refer to typical or av- 
erage values of the respective quantities 3 . Note that the 
width w of the peak in P(a) is, close to the jamming 
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Fig. 5: Illustration of Wyart 's argument [25] for u»/uj_: Left: 
a floppy patch of material, obtained by cutting bonds on the 
outer edge, in which all contacts have exactly a = n/2 upon 
distortion, in accord with the definition of a floppy mode [23]. 
Right: a weakly distorted floppy mode (also called trial soft 
mode) can be thought of as a floppy mode that is distorted 
elastically on a scale £*. Accordingly all angles a are slightly 
different from tt/2, as indicated in the figure. The dashed lines 
denote the relative displacement pairs of contacting particles, 
marked by the solid line connecting their centers. 



transition 

U\ 



ij/u±,ij if M||. 



/u_L, because 



tt/21 



roughly w ^ 

The question is now, whether the linear response follows 
this prediction for the soft modes, for our two families of 
networks. The insets of fig. 4 show that the scaling be- 
havior (6) is consistent with our data for the width w 
of the peak of P(a) for packing-derived networks under 
shear, and for randomly cut networks under either com- 
pression and shear deformations. The P{a) for compres- 
sion of packing-derived networks is the exceptional case. 
The peak of P(a) does not grow as much, and a substan- 
tial shoulder for large a remains even close to jamming: 
the tendency for particles to move towards each other re- 
mains much more prominent under compression. 

Fraction of Compressed Bonds. In order to clar- 
ify the significance of the large-a shoulder, let us analyze 
the fraction of significantly compressed bonds. Intuitively, 
it is clear that this fraction should be at the root of the 
difference between randomly cut networks, whose com- 
pression modulus K vanishes near jamming, and packing- 
derived networks whose K does not. Indeed, suppose we 
compress a packing-derived network. When a finite frac- 
tion of the bonds gets shortened with a finite fraction of 
the strain e, then K will be proportional to the bond 
strength k — this simply follows from the expression for 
the energy change AE, eq. (1). 

To quantify this, we define the fraction p comp of bonds 



3 In earlier work [7], we have argued that the scaling um/ux ~ Az 
can also be understood by balancing the first and second terms in 



the energy expansion (eq. (1)) which yields the scaling Un/ux ~ v 7 ^, 
with 5 is the typical overlap. For jammed packings, where the pre- 
stress term is taken into account, this result is consistent with (6) 
in view of the scaling Az ~ \J~5. However, as we show here, even if 
the pre-stress term is ignored in the dynamical matrix, very similar 
scaling is obtained. 
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Fig. 6: Scaling of (a) p C om P , (b) p s trict, and (c) py, as a function 
of z for compression of packing-derived networks (symbols) and 
randomly cut networks (curves), (d) The packing data on a 
log scale, emphasizing the rapid rise of p comp at small Az. The 
dotted lines mark exponents 0.5 and 0.65, to guide the eye — 
there is no clean scaling. 

whose local response has a > 37r/4, i.e., contact pairs 
which upon compression move more towards each other 
than they move sideways (un < — \u±\ < 0). If the PDF 
P(a) was governed by a single scale, the observed scaling 
of the width of the peak of P{a) ~ Az in accord with (6), 
suggests that /3 C omp ~ Az near jamming. 

In fig. 6 we compare the scaling of p comp for both our 
types of networks. The P(a) for random networks can be 
described by a single scale l/l* ~ w ~ Az (fig. 4b), and in- 
deed the corresponding p comp is linear in Az (fig. 6a). For 
packing-derived networks close to jamming, p comp rises 
more rapidly than linearly, and is much larger than for 
the randomly cut networks. This shows that under com- 
pression of packings a significant fraction of the contacts 
remains non-sliding and that single-parameter scaling does 
not apply — indeed, while our randomly cut networks are 
consistent with a linear variation of p comp , if we fit our 
data for packing-derived networks to a power law form 
Pcomp ~ (Az)^, we do not find a clear scaling (£ m 0.65, 
but only over 1 decade in Az). 

In principle, many of the bonds with a > 37r/4 could 
have anomalously small «|| — to check that this is not 
the case, we have also investigated pstrict , the fraction of 
bonds whose local response has a > 37r/4 and u\\ < —erij, 
and pu , the fraction of bonds whose local response has 
U|| < —erij. The latter condition can be phrased as the 
fraction of bonds that are more compressed than they 
would be if the response were affine. As shown in fig. 6b- 
c, these measures of compressed bond fractions are also 
much larger for compression of packing-derived networks. 

For compression of packing-derived networks, since K/k 
remains finite for Az — > 0, one should expect a finite frac- 
tion of bonds with it|| of order e — consistent with this p|| 



remains finite in this case. Although the P(a) in fig. 4a 
do not appear to be governed by a single scale, a tentative 
argument can be given why the rise in p com p is steeper 
than linear for small Az: Assume the typical u± is still 
of the order ej\J Az, as is the case for compression of ran- 
domly cut networks (from combining eq. (6) with eq. (1) 
and K ~ fcAz). Then, the relevant scale in P(a) would be 
set by uu/u± ~ \J Az, and one would find p comp ~ \J Az. 
As expected, we do not find such a clear scaling in fig. 6d, 
but the rapid initial rise is clearly visible. 

In conclusion, we find that the non-affine displacements 
in random spring networks and sheared jammed packings 
all share the same simple scalings of P(a), as well as hav- 
ing clastic moduli which scale as fcAz, where fc denotes the 
bond stiffness. The response of jammed packings to com- 
pression is the exceptional case: P(a) has more structure 
than a single peak, naive scaling breaks down and K ~ k. 

Interpretation in terms of the space of force net- 
works. — We finally briefly discuss these issues within 
the framework developed by Wyart [20,21] for the response 
of frictionless granular packings. For a network consisting 
of N particles and zTV/2 contacts, any imposed deforma- 
tion can be expressed in terms of the change of the rest 
lengths of some bonds in the network. After perturbing 
one or more bonds, for example in a way which corre- 
sponds to a global shear or compression of the packing, 
there will be an energy minimization involving the dN de- 
grees of freedom (displacements Uj). Hence, the space of 
responses to perturbations that cost energy has dimension 
zN/2 — dN — AzN/2. An equivalent way to view this is 
that after perturbing the rest lengths of the bonds, the 
particles will move so as to satisfy the dN local equations 
of force balance. Therefore the force response network can 
be expressed in a basis {f^} of the AzA/2-dimensional 
solution space T of of the force balance equations. 

The force space thus denned is very similar to the so- 
lution space of the force network ensemble [26-28], where 
one studies the space of allowed force configurations, T ine 
for a given contact geometry and externally imposed pres- 
sure. Let us define the extended force network ensemble, 
as the ensemble of all allowed force configurations, without 
the constraint that the pressure be fixed [28]. This force 
space is precisely the AzA/2-dimcnsional space spanned 
by the orthonormal basis {fW} defined above. 

Now, if we fix the pressure, this leads to an additional 
constraint. By a simple rotation in force space it is possi- 
ble to choose the {f^} such that f^ precisely gives the 
direction of increasing pressure, so that all other base vec- 
tors are perpendicular to the pressure direction [28] - 
the force ensemble with fixed pressure simply results from 
projecting out the f*- 1 -* direction from T . 

Suppose we externally impose changes in the rest 
lengths of the bonds, denoted by y — for a compression, 
we may for example increase all rest lengths. Wyart [20] 
then shows that the energy change corresponding to such 
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external forcing can be expressed as 
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part of this work was done, for its hospitality. 



(7) 



where (-|-) denotes the scalar product in force space. If we 
consider a deformation y of which the components are of 
order e then in general we may assume that there will be 
no correlation between y and the force space T . The domi- 
nant contribution to the squared inner product in eq. (7) is 



then Em(/> 



W\2„2 



X)m e 2 /-^ ~ e 2 , and summing these 



over all basis vectors gives AE ~ NAze 2 , making the 
energy extensive and proportional to the distance to the 
jamming transition times the square of the strain. This is 
the case for example if y represents a shear deformation, 
so that G ~ Az. 

In this scenario, however, the response to compression is 
an exceptional case. Compression amounts to the special 
situation that y is essentially pointing in the same direc- 
tion as the basis vector fW, which we chose to be in the 
direction of increasing pressure. In this case, all terms in 
the inner product Yl m f™ Vm are positive, and of order 
e/y/N, so that (f (1) |y) 2 ~ Ne 2 . Now the contribution of 
the first term in eq. (7) is already extensive by itself, and 
independent of Az. In this picture, the compression re- 
sponse is anomalous since it corresponds to an alignment 
with a special direction in the force space; in other words, 
from this point of view too, it is best to think of compres- 
sion, and not shear, as being the anomalous response! 

Conclusion. — Since the jamming transition is about 
loss of rigidity, it is natural to compare the elasticity of disc 
packings to rigidity percolation on a random network [8] . 
From this reference point, what is special about jammed 
packings is that they have an anomalously large resistance 
to compression. In addition, the response of packings to 
compression docs not follow the simple scaling relations 
that govern the particles' relative displacements under 
shear. This interpretation should be contrasted with the 
more commonly stated conclusion from comparing to ef- 
fective medium theory, which would be that the resistance 
to shear is anomalously small [2, 10]. 

The anomalous behavior of the bulk modulus has a 
purely geometric origin, as we have shown by numerically 
extracting it from the dynamical matrix, and interpret- 
ing the results in terms of the solution space of the force 
network ensemble, both of which are purely geometric: It 
is all a matter of how purely repulsive particles arrange 
themselves when they are pressed together. 
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